﻿#include <gtest/gtest.h>

#include <random>
#include <time.h>

#include <xy/geom/vec.h>
#include <xy/math/poly.h>

using namespace xy;
using namespace xy::math;
using namespace xy::geom;

TEST(TestPolynomial, Init)
{
    poly<double> p1({1, 2, 3});
    poly<double> p2({1, 2, 3, 4, 5});

    std::vector<double> coeffs = {1, 2, 3, 4, 5};
    poly<double> p3(std::move(coeffs));
    EXPECT_EQ(p2, p3);

    poly<double> p4(1);
    p4 = p1;
    EXPECT_EQ(p1, p4);

    poly<double> p5(0);
    p5 = std::move(p2);
    EXPECT_EQ(p3, p5);

    poly<double> p6 = std::move(p1);
    EXPECT_EQ(p4, p6);
}

TEST(TestPolynomial, Initvec3)
{
    poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
    poly<vec3, double> p2(
        {vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});

    std::vector<vec3> coeffs = {vec3{1.0, 1.0, 1.0}, vec3{2.0, 2.0, 2.0}, vec3{3.0, 3.0, 3.0}, vec3{4.0, 4.0, 4.0},
                                vec3{5.0, 5.0, 5.0}};
    poly<vec3, double> p3(std::move(coeffs));
    EXPECT_EQ(p2, p3);

    poly<vec3, double> p4(1);
    p4 = p1;
    EXPECT_EQ(p1, p4);

    poly<vec3, double> p5(0);
    p5 = std::move(p2);
    EXPECT_EQ(p3, p5);

    poly<vec3, double> p6 = std::move(p1);
    EXPECT_EQ(p4, p6);
}

TEST(TestPiecewise_Polynomial, Init)
{
    std::mt19937 gen(std::random_device{}());
    std::uniform_real_distribution<double> dis(-10, 10);

    {
        poly<double> p1({1, 2, 3});
        poly<double> p2 = p1;

        ppoly<double> pp1(std::move(p1));
        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            EXPECT_EQ(pp1.value(t), p2.value(t));
        }
    }

    {
        ppoly<double> pp1({1, 2, 3});
        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            EXPECT_EQ(pp1.value(t), 0);
        }
    }

    {
        poly<double> p1({1, 2, 3});
        poly<double> p2({1, 2, 3, 4, 5});
        poly<double> p3({1, 2, 3, 4, 5, 6});

        std::vector<double> breaks = {-3, 3};
        std::vector<poly<double>> pieces = {p1, p2, p3};
        ppoly<double> pp1(std::move(breaks), std::move(pieces));

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp1.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp1.value(t), p2.value(t));
            else
                EXPECT_EQ(pp1.value(t), p3.value(t));
        }

        ppoly<double> pp2 = pp1;
        ppoly<double> pp3 = std::move(pp1);

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp2.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp2.value(t), p2.value(t));
            else
                EXPECT_EQ(pp2.value(t), p3.value(t));
        }

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp3.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp3.value(t), p2.value(t));
            else
                EXPECT_EQ(pp3.value(t), p3.value(t));
        }

        ppoly<double> pp4({1, 2, 3});
        pp4 = pp2;
        ppoly<double> pp5({1, 2, 3});
        pp5 = std::move(pp2);
        EXPECT_EQ(pp4, pp5);
    }
}

TEST(TestPiecewise_Polynomial, Initvec3)
{
    std::mt19937 gen(std::random_device{}());
    std::uniform_real_distribution<double> dis(-10, 10);

    {
        poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> p2 = p1;

        ppoly<vec3, double> pp1(std::move(p1));
        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            EXPECT_EQ(pp1.value(t), p2.value(t));
        }
    }

    {
        ppoly<vec3, double> pp1({1.0, 2.0, 3.0});
        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            EXPECT_EQ(pp1.value(t), vec3(0.0, 0.0, 0.0));
        }
    }

    {
        poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> p2(
            {vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        poly<vec3, double> p3({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0),
                               vec3(5.0, 5.0, 5.0), vec3(6.0, 6.0, 6.0)});

        std::vector<double> breaks = {-3, 3};
        std::vector<poly<vec3, double>> pieces = {p1, p2, p3};
        ppoly<vec3, double> pp1(std::move(breaks), std::move(pieces));

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp1.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp1.value(t), p2.value(t));
            else
                EXPECT_EQ(pp1.value(t), p3.value(t));
        }

        ppoly<vec3, double> pp2 = pp1;
        ppoly<vec3, double> pp3 = std::move(pp1);

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp2.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp2.value(t), p2.value(t));
            else
                EXPECT_EQ(pp2.value(t), p3.value(t));
        }

        for (int i = 0; i < 1000; i++)
        {
            double t = dis(gen);
            if (t < -3)
                EXPECT_EQ(pp3.value(t), p1.value(t));
            else if (t < 3)
                EXPECT_EQ(pp3.value(t), p2.value(t));
            else
                EXPECT_EQ(pp3.value(t), p3.value(t));
        }

        ppoly<vec3, double> pp4({1.0, 2.0, 3.0});
        pp4 = pp2;
        ppoly<vec3, double> pp5({1.0, 2.0, 3.0});
        pp5 = std::move(pp2);
        EXPECT_EQ(pp4, pp5);
    }
}

TEST(TestPolynomial, Four)
{
    poly<double> p1({1, 2, 3});
    poly<double> p2({1, 2, 3, 4, 5});

    {
        poly<double> p3 = p1 + p2;
        poly<double> r3({2, 4, 6, 4, 5});
        EXPECT_EQ(p3, r3);

        poly<double> p4 = p1 - p2;
        poly<double> r4({0, 0, 0, -4, -5});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<double> p3 = p2 + p1;
        poly<double> r3({2, 4, 6, 4, 5});
        EXPECT_EQ(p3, r3);

        poly<double> p4 = p2 - p1;
        poly<double> r4({0, 0, 0, 4, 5});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<double> p3 = p2;
        p3 += p1;
        poly<double> r3({2, 4, 6, 4, 5});
        EXPECT_EQ(p3, r3);

        poly<double> p4 = p2;
        p4 -= p1;
        poly<double> r4({0, 0, 0, 4, 5});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<double> p5 = p1 * 2;
        poly<double> r5({2, 4, 6});
        EXPECT_EQ(p5, r5);

        poly<double> p6 = p1 / 2;
        poly<double> r6({0.5, 1, 1.5});
        EXPECT_EQ(p6, r6);
    }

    {
        poly<double> p5 = p1;
        p5 *= 2;
        poly<double> r5({2, 4, 6});
        EXPECT_EQ(p5, r5);

        poly<double> p6 = p1;
        p6 /= 2;
        poly<double> r6({0.5, 1, 1.5});
        EXPECT_EQ(p6, r6);
    }

    {
        poly<double> poly1({1, 2, 3});
        poly<double> poly2({2, 3, 4});
        poly<double> poly3({2, 7, 16, 17, 12});
        poly<double> poly4({1, 3, 2});
        poly<double> poly5({23.0 / 4, -0.5, 6});

        EXPECT_EQ(poly1 * poly2, poly3);
        EXPECT_EQ(poly1 *= poly2, poly3);
        EXPECT_EQ(poly3 / poly4, poly5);
        EXPECT_EQ(poly3 /= poly4, poly5);
    }
}

TEST(TestPolynomial, Fourvec3)
{
    poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
    poly<vec3, double> p2(
        {vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});

    {
        poly<vec3, double> p3 = p1 + p2;
        poly<vec3, double> r3(
            {vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        EXPECT_EQ(p3, r3);

        poly<vec3, double> p4 = p1 - p2;
        poly<vec3, double> r4({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(-4.0, -4.0, -4.0),
                               vec3(-5.0, -5.0, -5.0)});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<vec3, double> p3 = p2 + p1;
        poly<vec3, double> r3(
            {vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        EXPECT_EQ(p3, r3);

        poly<vec3, double> p4 = p2 - p1;
        poly<vec3, double> r4(
            {vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<vec3, double> p3 = p2;
        p3 += p1;
        poly<vec3, double> r3(
            {vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        EXPECT_EQ(p3, r3);

        poly<vec3, double> p4 = p2;
        p4 -= p1;
        poly<vec3, double> r4(
            {vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        EXPECT_EQ(p4, r4);
    }

    {
        poly<vec3, double> p5 = p1 * 2;
        poly<vec3, double> r5({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0)});
        EXPECT_EQ(p5, r5);

        poly<vec3, double> p6 = p1 / 2;
        poly<vec3, double> r6({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5)});
        EXPECT_EQ(p6, r6);
    }

    {
        poly<vec3, double> p5 = p1;
        p5 *= 2;
        poly<vec3, double> r5({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0)});
        EXPECT_EQ(p5, r5);

        poly<vec3, double> p6 = p1;
        p6 /= 2;
        poly<vec3, double> r6({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5)});
        EXPECT_EQ(p6, r6);
    }

    {
        poly<vec3, double> poly1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> poly2({vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0)});
        poly<vec3, double> poly3({vec3(2.0, 2.0, 2.0), vec3(7.0, 7.0, 7.0), vec3(16.0, 16.0, 16.0),
                                  vec3(17.0, 17.0, 17.0), vec3(12.0, 12.0, 12.0)});
        poly<vec3, double> poly4({vec3(1.0, 1.0, 1.0), vec3(3.0, 3.0, 3.0), vec3(2.0, 2.0, 2.0)});
        poly<vec3, double> poly5({vec3(23.0 / 4, 23.0 / 4, 23.0 / 4), vec3(-0.5, -0.5, -0.5), vec3(6.0, 6.0, 6.0)});

        EXPECT_EQ(poly1 * poly2, poly3);
        EXPECT_EQ(poly1 *= poly2, poly3);
        EXPECT_EQ(poly3 / poly4, poly5);
        EXPECT_EQ(poly3 /= poly4, poly5);
    }
}

TEST(TestPiecewise_Polynomial, Four)
{
    poly<double> p1({1, 2, 3});
    poly<double> p2({1, 2, 3, 4, 5});
    ppoly<double> pp1({3}, {p1, p2});
    ppoly<double> pp2({3}, {p2, p1});

    {
        ppoly<double> pp3 = pp1 + pp2;
        ppoly<double> r3({3}, {poly<double>({2, 4, 6, 4, 5}), poly<double>({2, 4, 6, 4, 5})});
        EXPECT_EQ(pp3, r3);

        ppoly<double> pp4 = pp1 - pp2;
        ppoly<double> r4({3}, {poly<double>({0, 0, 0, -4, -5}), poly<double>({0, 0, 0, 4, 5})});
        EXPECT_EQ(pp4, r4);
    }

    {
        ppoly<double> pp3 = pp2 + pp1;
        ppoly<double> r3({3}, {poly<double>({2, 4, 6, 4, 5}), poly<double>({2, 4, 6, 4, 5})});
        EXPECT_EQ(pp3, r3);

        ppoly<double> pp4 = pp2 - pp1;
        ppoly<double> r4({3}, {poly<double>({0, 0, 0, 4, 5}), poly<double>({0, 0, 0, -4, -5})});
        EXPECT_EQ(pp4, r4);
    }
    {
        ppoly<double> pp3 = pp2;
        pp3 += pp1;
        ppoly<double> r3({3}, {poly<double>({2, 4, 6, 4, 5}), poly<double>({2, 4, 6, 4, 5})});
        EXPECT_EQ(pp3, r3);

        ppoly<double> pp4 = pp2;
        pp4 -= pp1;
        ppoly<double> r4({3}, {poly<double>({0, 0, 0, 4, 5}), poly<double>({0, 0, 0, -4, -5})});
        EXPECT_EQ(pp4, r4);
    }

    {
        ppoly<double> pp5 = pp1 * 2;
        ppoly<double> r5({3}, {poly<double>({2, 4, 6}), poly<double>({2, 4, 6, 8, 10})});
        EXPECT_EQ(pp5, r5);

        ppoly<double> pp6 = pp1 / 2;
        ppoly<double> r6({3}, {poly<double>({0.5, 1, 1.5}), poly<double>({0.5, 1, 1.5, 2, 2.5})});
        EXPECT_EQ(pp6, r6);
    }

    {
        ppoly<double> pp5 = pp1;
        pp5 *= 2;
        ppoly<double> r5({3}, {poly<double>({2, 4, 6}), poly<double>({2, 4, 6, 8, 10})});
        EXPECT_EQ(pp5, r5);

        ppoly<double> pp6 = pp1;
        pp6 /= 2;
        ppoly<double> r6({3}, {poly<double>({0.5, 1, 1.5}), poly<double>({0.5, 1, 1.5, 2, 2.5})});
        EXPECT_EQ(pp6, r6);
    }

    {
        poly<double> poly1({1, 2, 3});
        poly<double> poly2({2, 3, 4});
        poly<double> poly3({2, 7, 16, 17, 12});
        poly<double> poly4({1, 3, 2});
        poly<double> poly5({23.0 / 4, -0.5, 6});

        ppoly<double> pp1({3}, {poly1, poly2});
        ppoly<double> pp2({3}, {poly2, poly1});
        ppoly<double> pp3({3}, {poly3, poly3});
        ppoly<double> pp4({3}, {poly4, poly4});
        ppoly<double> pp5({3}, {poly5, poly5});

        EXPECT_EQ(pp1 * pp2, pp3);
        EXPECT_EQ(pp1 *= pp2, pp3);
        EXPECT_EQ(pp3 / pp4, pp5);
        EXPECT_EQ(pp3 /= pp4, pp5);
    }
}

TEST(TestPiecewise_Polynomial, Fourvec3)
{
    poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
    poly<vec3, double> p2(
        {vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
    ppoly<vec3, double> pp1({3}, {p1, p2});
    ppoly<vec3, double> pp2({3}, {p2, p1});

    {
        ppoly<vec3, double> pp3 = pp1 + pp2;
        ppoly<vec3, double> r3({3}, {poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)}),
                                     poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)})});
        EXPECT_EQ(pp3, r3);

        ppoly<vec3, double> pp4 = pp1 - pp2;
        ppoly<vec3, double> r4({3}, {poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(-4.0, -4.0, -4.0), vec3(-5.0, -5.0, -5.0)}),
                                     poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)})});
        EXPECT_EQ(pp4, r4);
    }

    {
        ppoly<vec3, double> pp3 = pp2 + pp1;
        ppoly<vec3, double> r3({3}, {poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)}),
                                     poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)})});
        EXPECT_EQ(pp3, r3);

        ppoly<vec3, double> pp4 = pp2 - pp1;
        ppoly<vec3, double> r4({3}, {poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)}),
                                     poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(-4.0, -4.0, -4.0), vec3(-5.0, -5.0, -5.0)})});
        EXPECT_EQ(pp4, r4);
    }
    {
        ppoly<vec3, double> pp3 = pp2;
        pp3 += pp1;
        ppoly<vec3, double> r3({3}, {poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)}),
                                     poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)})});
        EXPECT_EQ(pp3, r3);
        ppoly<vec3, double> pp4 = pp2;
        pp4 -= pp1;
        ppoly<vec3, double> r4({3}, {poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)}),
                                     poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0),
                                                         vec3(-4.0, -4.0, -4.0), vec3(-5.0, -5.0, -5.0)})});
        EXPECT_EQ(pp4, r4);
    }

    {
        ppoly<vec3, double> pp5 = pp1 * 2;
        ppoly<vec3, double> r5({3},
                               {poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0)}),
                                poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                    vec3(8.0, 8.0, 8.0), vec3(10.0, 10.0, 10.0)})});
        EXPECT_EQ(pp5, r5);

        ppoly<vec3, double> pp6 = pp1 / 2;
        ppoly<vec3, double> r6({3},
                               {poly<vec3, double>({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5)}),
                                poly<vec3, double>({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5),
                                                    vec3(2.0, 2.0, 2.0), vec3(2.5, 2.5, 2.5)})});
        EXPECT_EQ(pp6, r6);
    }

    {
        ppoly<vec3, double> pp5 = pp1;
        pp5 *= 2;
        ppoly<vec3, double> r5({3},
                               {poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0)}),
                                poly<vec3, double>({vec3(2.0, 2.0, 2.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0),
                                                    vec3(8.0, 8.0, 8.0), vec3(10.0, 10.0, 10.0)})});
        EXPECT_EQ(pp5, r5);

        ppoly<vec3, double> pp6 = pp1;
        pp6 /= 2;
        ppoly<vec3, double> r6({3},
                               {poly<vec3, double>({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5)}),
                                poly<vec3, double>({vec3(0.5, 0.5, 0.5), vec3(1.0, 1.0, 1.0), vec3(1.5, 1.5, 1.5),
                                                    vec3(2.0, 2.0, 2.0), vec3(2.5, 2.5, 2.5)})});
        EXPECT_EQ(pp6, r6);
    }

    {
        poly<vec3, double> poly1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> poly2({vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0)});
        poly<vec3, double> poly3({vec3(2.0, 2.0, 2.0), vec3(7.0, 7.0, 7.0), vec3(16.0, 16.0, 16.0),
                                  vec3(17.0, 17.0, 17.0), vec3(12.0, 12.0, 12.0)});
        poly<vec3, double> poly4({vec3(1.0, 1.0, 1.0), vec3(3.0, 3.0, 3.0), vec3(2.0, 2.0, 2.0)});
        poly<vec3, double> poly5({vec3(23.0 / 4, 23.0 / 4, 23.0 / 4), vec3(-0.5, -0.5, -0.5), vec3(6.0, 6.0, 6.0)});

        ppoly<vec3, double> pp1({3}, {poly1, poly2});
        ppoly<vec3, double> pp2({3}, {poly2, poly1});
        ppoly<vec3, double> pp3({3}, {poly3, poly3});
        ppoly<vec3, double> pp4({3}, {poly4, poly4});
        ppoly<vec3, double> pp5({3}, {poly5, poly5});

        EXPECT_EQ(pp1 * pp2, pp3);
        EXPECT_EQ(pp1 *= pp2, pp3);
        EXPECT_EQ(pp3 / pp4, pp5);
        EXPECT_EQ(pp3 /= pp4, pp5);
    }
}

TEST(TestPiecewise_Polynomial, Op)
{
    // 推入节点并设置多项式
    {
        poly<double> p1({1, 2, 3});
        poly<double> p2({1, 2, 3, 4, 5});
        poly<double> p3({2, 7, 16, 17, 12});
        poly<double> p4({1, 3, 2});
        ppoly<double> pp1({3}, {p1, p2});
        ppoly<double> pp2({1, 3, 5}, {p4, p1, p2, p3});

        pp1.push_back(5) = p3;
        pp1.push_front(1) = p4;

        EXPECT_EQ(pp1, pp2);
    }

    // 求导，积分
    {
        poly<double> p1({1, 2, 3});
        poly<double> p2({1, 2, 3, 4, 5});
        poly<double> p3({2, 7, 16, 17, 12});
        poly<double> p4({1, 3, 2});
        ppoly<double> pp1({-3, 0, 3}, {p1, p2, p3, p4});
        ppoly<double> pp2({-3, 0, 3}, {p1.derivative(), p2.derivative(), p3.derivative(), p4.derivative()});
        ppoly<double> pp3({-3, 0, 3}, {p1.integrate(), p2.integrate(), p3.integrate(), p4.integrate()});

        EXPECT_EQ(pp1.derivative(), pp2);
        EXPECT_EQ(pp1.integrate(), pp3);

        EXPECT_DOUBLE_EQ(pp1.integrate(-5, -4), p1.integrate(-5, -4));
        EXPECT_DOUBLE_EQ(pp1.integrate(-5, -1), p1.integrate(-5, -3) + p2.integrate(-3, -1));
        EXPECT_DOUBLE_EQ(pp1.integrate(-5, 2), p1.integrate(-5, -3) + p2.integrate(-3, 0) + p3.integrate(0, 2));
        EXPECT_DOUBLE_EQ(pp1.integrate(-5, 5),
                         p1.integrate(-5, -3) + p2.integrate(-3, 0) + p3.integrate(0, 3) + p4.integrate(3, 5));

        EXPECT_DOUBLE_EQ(pp1.integrate(-2, -1), p2.integrate(-2, -1));
        EXPECT_DOUBLE_EQ(pp1.integrate(-2, 2), p2.integrate(-2, 0) + p3.integrate(0, 2));
        EXPECT_DOUBLE_EQ(pp1.integrate(-2, 5), p2.integrate(-2, 0) + p3.integrate(0, 3) + p4.integrate(3, 5));

        EXPECT_DOUBLE_EQ(pp1.integrate(1, 2), p3.integrate(1, 2));
        EXPECT_DOUBLE_EQ(pp1.integrate(1, 5), p3.integrate(1, 3) + p4.integrate(3, 5));

        EXPECT_DOUBLE_EQ(pp1.integrate(3, 5), p4.integrate(3, 5));
    }

    // 连续性测试
    {
        poly<double> p1({1, 1, 1, 1});
        poly<double> p2({-2, 2, 2, 2});
        poly<double> p3({18, -14, 3, 3});
        poly<double> p4({-9, 13, -6, 4});
        poly<double> p5({147, 6, -18, 5});
        poly<double> p6({947, -289, 4, 6});

        ppoly<double> pp1({1, 2, 3, 4, 5}, {p1, p2, p3, p4, p5, p6});
        std::vector<int> continuity = {0, 1, 2, 0, 1};
        for (int i = 0; i < 5; i++)
            EXPECT_EQ(pp1.continuity(i), continuity[i]);

        EXPECT_EQ(pp1.continuity(), 0);
    }

    // 合并测试
    {
        std::mt19937 generator(std::random_device{}());
        std::uniform_real_distribution<double> uniform(0, 1.0);
        std::uniform_int_distribution<int> int_uniform(1, 10);

        auto f = [&]() {
            int d = int_uniform(generator) % 5;
            poly<double> p(d);
            for (int i = 0; i < d + 1; i++)
                p.coefficient(i) = uniform(generator);
            return p;
        };

        std::vector<double> x(int_uniform(generator));
        std::vector<double> y(int_uniform(generator));

        x[0] = uniform(generator);
        for (int i = 1; i < x.size(); i++)
            x[i] = x[i - 1] + uniform(generator);

        y[0] = uniform(generator);
        for (int i = 1; i < y.size(); i++)
            y[i] = y[i - 1] + uniform(generator);

        ppoly<double> pp1(std::move(x));
        ppoly<double> pp2(std::move(y));

        for (int i = 0; i < pp1.poly_num(); i++)
            pp1[i] = f();

        for (int i = 0; i < pp2.poly_num(); i++)
            pp2[i] = f();

        auto pp3 = merge(pp1, pp2);

        double a = pp3.intervals().front() - 1;
        double b = pp3.intervals().back() + 1;
        double h = (b - a) / 1000;

        for (double x = a; x <= b; x += h)
        {
            double y1 = pp1.value(x);
            double y2 = pp2.value(x);
            double y3 = pp3.value(x);
            EXPECT_FLOAT_EQ(y3, y1 + y2);
        }
    }
}

TEST(TestPiecewise_Polynomial, Opvec3)
{
    // 推入节点并设置多项式
    {
        poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> p2({vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0)});
        poly<vec3, double> p3({vec3(2.0, 2.0, 2.0), vec3(7.0, 7.0, 7.0), vec3(16.0, 16.0, 16.0), vec3(17.0, 17.0, 17.0),
                               vec3(12.0, 12.0, 12.0)});
        poly<vec3, double> p4({vec3(1.0, 1.0, 1.0), vec3(3.0, 3.0, 3.0), vec3(2.0, 2.0, 2.0)});
        ppoly<vec3, double> pp1({3}, {p1, p2});
        ppoly<vec3, double> pp2({1, 3, 5}, {p4, p1, p2, p3});

        pp1.push_back(5) = p3;
        pp1.push_front(1) = p4;

        EXPECT_EQ(pp1, pp2);
    }

    // 求导，积分
    {
        poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> p2(
            {vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0), vec3(5.0, 5.0, 5.0)});
        poly<vec3, double> p3({vec3(2.0, 2.0, 2.0), vec3(7.0, 7.0, 7.0), vec3(16.0, 16.0, 16.0), vec3(17.0, 17.0, 17.0),
                               vec3(12.0, 12.0, 12.0)});
        poly<vec3, double> p4({vec3(1.0, 1.0, 1.0), vec3(3.0, 3.0, 3.0), vec3(2.0, 2.0, 2.0)});
        ppoly<vec3, double> pp1({-3, 0, 3}, {p1, p2, p3, p4});
        ppoly<vec3, double> pp2({-3, 0, 3}, {p1.derivative(), p2.derivative(), p3.derivative(), p4.derivative()});
        ppoly<vec3, double> pp3({-3, 0, 3}, {p1.integrate(), p2.integrate(), p3.integrate(), p4.integrate()});

        EXPECT_EQ(pp1.derivative(), pp2);
        EXPECT_EQ(pp1.integrate(), pp3);

        EXPECT_EQ(pp1.integrate(-5, -4), p1.integrate(-5, -4));
        EXPECT_EQ(pp1.integrate(-5, -1), p1.integrate(-5, -3) + p2.integrate(-3, -1));
        EXPECT_EQ(pp1.integrate(-5, 2), p1.integrate(-5, -3) + p2.integrate(-3, 0) + p3.integrate(0, 2));
        EXPECT_EQ(pp1.integrate(-5, 5),
                  p1.integrate(-5, -3) + p2.integrate(-3, 0) + p3.integrate(0, 3) + p4.integrate(3, 5));

        EXPECT_EQ(pp1.integrate(-2, -1), p2.integrate(-2, -1));
        EXPECT_EQ(pp1.integrate(-2, 2), p2.integrate(-2, 0) + p3.integrate(0, 2));
        EXPECT_EQ(pp1.integrate(-2, 5), p2.integrate(-2, 0) + p3.integrate(0, 3) + p4.integrate(3, 5));

        EXPECT_EQ(pp1.integrate(1, 2), p3.integrate(1, 2));
        EXPECT_EQ(pp1.integrate(1, 5), p3.integrate(1, 3) + p4.integrate(3, 5));

        EXPECT_EQ(pp1.integrate(3, 5), p4.integrate(3, 5));
    }

    // 连续性测试
    {
        poly<vec3, double> p1({vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0)});
        poly<vec3, double> p2({vec3(-2.0, -2.0, -2.0), vec3(2.0, 2.0, 2.0), vec3(2.0, 2.0, 2.0), vec3(2.0, 2.0, 2.0)});
        poly<vec3, double> p3(
            {vec3(18.0, 18.0, 18.0), vec3(-14.0, -14.0, -14.0), vec3(3.0, 3.0, 3.0), vec3(3.0, 3.0, 3.0)});
        poly<vec3, double> p4(
            {vec3(-9.0, -9.0, -9.0), vec3(13.0, 13.0, 13.0), vec3(-6.0, -6.0, -6.0), vec3(4.0, 4.0, 4.0)});
        poly<vec3, double> p5(
            {vec3(147.0, 147.0, 147.0), vec3(6.0, 6.0, 6.0), vec3(-18.0, -18.0, -18.0), vec3(5.0, 5.0, 5.0)});
        poly<vec3, double> p6(
            {vec3(947.0, 947.0, 947.0), vec3(-289.0, -289.0, -289.0), vec3(4.0, 4.0, 4.0), vec3(6.0, 6.0, 6.0)});

        ppoly<vec3, double> pp1({1, 2, 3, 4, 5}, {p1, p2, p3, p4, p5, p6});
        std::vector<int> continuity = {0, 1, 2, 0, 1};
        for (int i = 0; i < 5; i++)
            EXPECT_EQ(pp1.continuity(i), continuity[i]);

        EXPECT_EQ(pp1.continuity(), 0);
    }

    // 合并测试
    {
        std::mt19937 generator(std::random_device{}());
        std::uniform_real_distribution<double> uniform(0, 1.0);
        std::uniform_int_distribution<int> int_uniform(1, 10);

        auto f = [&]() {
            int d = int_uniform(generator) % 5;
            poly<vec3, double> p(d);
            for (int i = 0; i < d + 1; i++)
                p.coefficient(i) = vec3(uniform(generator), uniform(generator), uniform(generator));
            return p;
        };

        std::vector<double> x(int_uniform(generator));
        std::vector<double> y(int_uniform(generator));

        x[0] = uniform(generator);
        for (int i = 1; i < x.size(); i++)
            x[i] = x[i - 1] + uniform(generator);

        y[0] = uniform(generator);
        for (int i = 1; i < y.size(); i++)
            y[i] = y[i - 1] + uniform(generator);

        ppoly<vec3, double> pp1(std::move(x));
        ppoly<vec3, double> pp2(std::move(y));

        for (int i = 0; i < pp1.poly_num(); i++)
            pp1[i] = f();

        for (int i = 0; i < pp2.poly_num(); i++)
            pp2[i] = f();

        auto pp3 = merge(pp1, pp2);

        double a = pp3.intervals().front() - 1;
        double b = pp3.intervals().back() + 1;
        double h = (b - a) / 1000;

        for (double x = a; x <= b; x += h)
        {
            vec3 y1 = pp1.value(x);
            vec3 y2 = pp2.value(x);
            vec3 y3 = pp3.value(x);
            EXPECT_EQ(y3, y1 + y2);
        }
    }
}

TEST(TestPolynomial, Evaluate)
{
    poly<double> poly1({1, 2, 3, 4});
    poly<double> poly2({2, 6, 12});
    poly<double> poly3({0, 1, 1, 1, 1});

    EXPECT_DOUBLE_EQ(poly1.value(2), 49);
    EXPECT_EQ(poly1.derivative(), poly2);
    EXPECT_EQ(poly1.integrate(), poly3);
    EXPECT_DOUBLE_EQ(poly1.integrate(1, 2), 26);
}

TEST(TestPolynomial, Evaluatevec3)
{
    poly<vec3, double> poly1({vec3(1.0, 1.0, 1.0), vec3(2.0, 2.0, 2.0), vec3(3.0, 3.0, 3.0), vec3(4.0, 4.0, 4.0)});
    poly<vec3, double> poly2({vec3(2.0, 2.0, 2.0), vec3(6.0, 6.0, 6.0), vec3(12.0, 12.0, 12.0)});
    poly<vec3, double> poly3(
        {vec3(0.0, 0.0, 0.0), vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0), vec3(1.0, 1.0, 1.0)});

    EXPECT_EQ(poly1.value(2), vec3(49.0, 49.0, 49.0));
    EXPECT_EQ(poly1.derivative(), poly2);
    EXPECT_EQ(poly1.integrate(), poly3);
    EXPECT_EQ(poly1.integrate(1, 2), vec3(26.0, 26.0, 26.0));
}

// 测试多项式算法
TEST(TestPolynomial, Algorithm)
{
    // 节点多项式
    {
        std::vector<double> c = {1, 2, 3};
        auto p = poly<double>::nodes(c);
        poly<double> p2({-6, 11, -6, 1});
        EXPECT_EQ(p, p2);
    }

    // Lagrange
    {
        std::vector<double> c = {1, 4, 2, 3};
        auto p = poly<double>::lagrange(c, 1);
        poly<double> p2({-6, 11, -6, 1});
        EXPECT_EQ(p, (p2 / 6));
    }

    // Legendre
    {
        auto p1 = poly<double>::legendre(1);
        auto p2 = poly<double>::legendre(2);
        auto p3 = poly<double>::legendre(3);

        // dbg(p1, p2, p3);

        EXPECT_EQ(p1, poly<double>({0, 1}));
        EXPECT_EQ(p2, poly<double>({-0.5, 0, 1.5}));
        EXPECT_EQ(p3, poly<double>({0, -1.5, 0, 2.5}));
    }

    // MonicLegendre
    {
        auto p1 = poly<double>::monic_legendre(1);
        auto p2 = poly<double>::monic_legendre(2);
        auto p3 = poly<double>::monic_legendre(3);

        // dbg(p1, p2, p3);

        EXPECT_EQ(p1, poly<double>({-0.5, 1.0}));
        EXPECT_EQ(p2, poly<double>({1.0 / 6, -1.0, 1.0}));
        EXPECT_EQ(p3, poly<double>({-0.05, 0.6, -1.5, 1.0}));
    }
}

// 测试多项式算法
TEST(TestPolynomial, Algorithmvec3)
{
    // 节点多项式
    {
        std::vector<double> c = {1, 2, 3};
        auto p = poly<vec3, double>::nodes(c);
        poly<vec3, double> p2(
            {vec3(-6.0, -6.0, -6.0), vec3(11.0, 11.0, 11.0), vec3(-6.0, -6.0, -6.0), vec3(1.0, 1.0, 1.0)});
        EXPECT_EQ(p, p2);
    }

    // Lagrange
    {
        std::vector<double> c = {1, 4, 2, 3};
        auto p = poly<vec3, double>::lagrange(c, 1);
        poly<vec3, double> p2(
            {vec3(-6.0, -6.0, -6.0), vec3(11.0, 11.0, 11.0), vec3(-6.0, -6.0, -6.0), vec3(1.0, 1.0, 1.0)});
        EXPECT_EQ(p, (p2 / 6));
    }

    // Legendre
    {
        auto p1 = poly<vec3, double>::legendre(1);
        auto p2 = poly<vec3, double>::legendre(2);
        auto p3 = poly<vec3, double>::legendre(3);

        // dbg(p1, p2, p3);

        EXPECT_EQ(p1, (poly<vec3, double>({vec3(0.0, 0.0, 0.0), vec3(1.0, 1.0, 1.0)})));
        EXPECT_EQ(p2, (poly<vec3, double>({vec3(-0.5, -0.5, -0.5), vec3(0.0, 0.0, 0.0), vec3(1.5, 1.5, 1.5)})));
        EXPECT_EQ(p3, (poly<vec3, double>(
                          {vec3(0.0, 0.0, 0.0), vec3(-1.5, -1.5, -1.5), vec3(0.0, 0.0, 0.0), vec3(2.5, 2.5, 2.5)})));
    }

    // MonicLegendre
    {
        auto p1 = poly<vec3, double>::monic_legendre(1);
        auto p2 = poly<vec3, double>::monic_legendre(2);
        auto p3 = poly<vec3, double>::monic_legendre(3);

        // dbg(p1, p2, p3);

        EXPECT_EQ(p1, (poly<vec3, double>({vec3(-0.5, -0.5, -0.5), vec3(1.0, 1.0, 1.0)})));
        EXPECT_EQ(p2,
                  (poly<vec3, double>({vec3(1.0 / 6, 1.0 / 6, 1.0 / 6), vec3(-1.0, -1.0, -1.0), vec3(1.0, 1.0, 1.0)})));
        EXPECT_EQ(p3, (poly<vec3, double>({vec3(-0.05, -0.05, -0.05), vec3(0.6, 0.6, 0.6), vec3(-1.5, -1.5, -1.5),
                                           vec3(1.0, 1.0, 1.0)})));
    }
}